import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mplHomework 2
Submission Deadline: 14.05.2024 - 23:59
- Your homework solution has to be handed in as a group solution via Moodle.
- Clearly state name and matriculation number of each student
1 Velocity Decomposition
Recall the decomposition of a velocity field \(\mathbf{v} \left( \mathbf{x} + \mathbf{r} \right)\) into translation, rotation and internal deformation:
\[ \mathbf{v} \left( \mathbf{x} + \mathbf{r} \right) = \mathbf{v} \left( \mathbf{x}\right) + \mathbf{w} \times \mathbf{r} + \mathbf{D} \cdot \mathbf{r} \]
Consider the following flow fields that describe plane motions (2-D):
A) \(\mathbf{v} \left( \mathbf{x} \right) = \begin{pmatrix}1 \\ 1 \\ 0 \end{pmatrix}\), B) \(\mathbf{v} \left( \mathbf{x} \right) = \frac{1}{x^2+y^2} \begin{pmatrix}x \\ y \\ 0 \end{pmatrix}\), C) \(\mathbf{v} \left( \mathbf{x} \right) = \begin{pmatrix}-y \\ x \\ 0 \end{pmatrix}\) and
D) \(\mathbf{v} \left( \mathbf{x} \right) = \frac{1}{x^2+y^2} \begin{pmatrix}-y \\ x \\ 0 \end{pmatrix}\).
Tasks
Task 1
Decompose the given velocity fields and calculate for each field: \(\mathbf{w}\) and \(\mathbf{D}\)
Task 2
Describe the motion produced by each velocity field.
Task 3
Sketch each flow field within a domain given by \(x,y \in [-4, 4]\times [-4, 4]\)
Task 4
Within the flow field, sketch how a volume element having vertices \((1, 1), (1, 2), (2, 2), (2,1)\) deforms.
2 Velocity Gradient, Divergence, Curl
The gradient, divergence and curl(or vorticity) of a two dimensional velocity field \(\mathbf{v} = \begin{pmatrix} v_x \\ v_y \end{pmatrix} \mathbf{e}_i\) are given as
\[ \text{grad}(\mathbf{v}) = \nabla \mathbf{v} = \begin{bmatrix} \partial_x v_x & \partial_y v_x \\ \partial_x v_y & \partial_y v_y \\ \end{bmatrix} \mathbf{e}_i \otimes \mathbf{e}_j \]
\[ \text{div}(\mathbf{v}) = \nabla \cdot \mathbf{v} = \partial_x v_x + \partial_y v_y \]
\[ \text{curl}(\mathbf{v}) = \text{vorticity} = \mathbf{\omega} = \nabla \times \mathbf{v} = (\partial_y v_x - \partial_x v_y) \mathbf{e}_3 \]
The goal of this problem is to compute and visulize the velocity field and its gradient, divergence and curl.
For the rest of the problem, consider the velocity field
\[ \mathbf{v} = \begin{pmatrix} \log(1+y^2) +\sin(x+y) +0.5(x-2.5) -0.7(x+2.5)\\ \log(1+x^2) -\cos(x-y) +0.5(y-2.5) -0.7(y+2.5) \end{pmatrix} \mathbf{e}_i \]
Setup
For the visualization, we shall use the python libraries numpy and matplotlib. Here, we first import the libraries
Let us define a grid of points within the domain \(x, y \in [-5, 5] \times [-5, 5]\)
x, y = np.meshgrid(np.linspace(-5,5,21),np.linspace(-5,5,21))Task1 - Velocity Field
Computation
Implement the computation for the components of the velocity vector at every point \((x, y)\) on the grid using numpy math functions.
#### YOUR SOLUTION HERE
v_x = np.zeros_like(x)
v_y = np.zeros_like(y)
####Visualization
We can now visualize the velocity field using a quiver plot
Code
plt.quiver(x,y,v_x, v_y)
plt.gca().set_aspect('equal')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Velocity field");
plt.xticks([-5, 0, 5]);
plt.yticks([-5, 0, 5]);
Task 2 - Velocity Gradient
The velocity gradient \(\nabla \mathbf{v}\) is a second order tensor, and hence cannot be visualized on a 2D plot. Here, we shall look at the individual components of the velocity gradient tensor.
Computation
Calculate the components of the velocity gradient tensor and implement the computation in the following
#### YOUR SOLUTION HERE
dvx_dx = np.zeros_like(x)
dvx_dy = np.zeros_like(y)
dvy_dx = np.zeros_like(x)
dvy_dy = np.zeros_like(y)
####Visualization
We can now visualize the components of the velocity gradient individually.
Code
figure, axs = plt.subplots(nrows=2, ncols=2, figsize=(9, 7))
min_val = np.min([np.min(dvx_dx), np.min(dvx_dy), np.min(dvy_dx), np.min(dvy_dy)])
max_val = np.max([np.max(dvx_dx), np.max(dvx_dy), np.max(dvy_dx), np.max(dvy_dy)])
axs[0, 0].imshow(dvx_dx, extent=(-5.5, 5.5, -5.5, 5.5), vmin=min_val, vmax=max_val, cmap='coolwarm', interpolation='spline16')
axs[0, 0].quiver(x,y,v_x, v_y)
axs[0, 0].set_xticks([-5, 0, 5])
axs[0, 0].set_yticks([-5, 0, 5])
axs[0, 0].set_xlabel("x")
axs[0, 0].set_ylabel("y")
axs[0, 0].set_title(r'$\partial_x v_x$')
axs[0, 1].imshow(dvx_dy, extent=(-5.5, 5.5, -5.5, 5.5), vmin=min_val, vmax=max_val, cmap='coolwarm', interpolation='spline16')
axs[0, 1].quiver(x,y,v_x, v_y)
axs[0, 1].set_xticks([-5, 0, 5])
axs[0, 1].set_yticks([-5, 0, 5])
axs[0, 1].set_xlabel("x")
axs[0, 1].set_ylabel("y")
axs[0, 1].set_title(r'$\partial_y v_x$')
axs[1, 0].imshow(dvy_dx, extent=(-5.5, 5.5, -5.5, 5.5), vmin=min_val, vmax=max_val, cmap='coolwarm', interpolation='spline16')
axs[1, 0].quiver(x,y,v_x, v_y)
axs[1, 0].set_xticks([-5, 0, 5])
axs[1, 0].set_yticks([-5, 0, 5])
axs[1, 0].set_xlabel("x")
axs[1, 0].set_ylabel("y")
axs[1, 0].set_title(r'$\partial_x v_y$')
im = axs[1, 1].imshow(dvy_dy, extent=(-5.5, 5.5, -5.5, 5.5), vmin=min_val, vmax=max_val, cmap='coolwarm', interpolation='spline16')
axs[1, 1].quiver(x,y,v_x, v_y)
axs[1, 1].set_xticks([-5, 0, 5])
axs[1, 1].set_yticks([-5, 0, 5])
axs[1, 1].set_xlabel("x");
axs[1, 1].set_ylabel("y");
axs[1, 1].set_title(r'$\partial_y v_y$');
figure.subplots_adjust(right=0.8)
cax,kw = mpl.colorbar.make_axes([ax for ax in axs.flat])
plt.colorbar(im, cax=cax, **kw)
plt.show()
Task 3 - Divergence of Velocity
Computation
Calculate the divergence of the velocity field using already computed quatities, and implement it in the following.
#### YOUR SOLUTION HERE
div_v = np.zeros_like(dvx_dx)
####Visualization
Plot the divergence of the velocity field and interpret the result
#### YOUR SOLUTION HERE
####Task 4 - Vorticity
Computation
Calculate the vorticity (or curl) of the velocity field using previously computed quantities and implement it in the following
#### YOUR SOLUTION HERE
curl_v = np.zeros_like(dvx_dx)
####Visualization
Plot the divergence of the velocity field and interpret the result.
#### YOUR SOLUTION HERE
####3 Venturi Effect
The figure shows a Venturi tube, which can be used to measure fluid flow. It comprises a cylindrical inlet with cross section \(A_1\) followed by a convergent entrance into a cylindrical throat with cross section \(A_2\) and a divergent outlet. The velocity of the fluid at inlet and the throat is \(\mathbf{v}_1\) and \(\mathbf{v}_2\) respectively. The flowing fluid has density \(\rho_a\) and the liquid in the U-shaped tube has a density \(\rho_f\).
Tasks
Task 1
State the Bernoulli equation and explain to what fluid flow scenarios can it be applied.
Task 2
Calculate the velocity of the flowing fluid in terms of the difference in height of the fluid \(\Delta h\) in the U-shaped tube and other given quantities given in the figure.
Task 3
If the velocity of the inflowing fluid is doubled, how would \(\Delta h\) change?
4 Spinning Fluid
Consider a stationary cylindrical glass with inner radius \(R\) and partially filled with water upto height \(H\). The ambient pressure is given by \(p_a\). The glass is now spun about its central axis with an angular velocity \(\Omega\), which causes the top surface of the water to form a parabolloid. A selected cross section of the stationary and spinning glass with fluid is shown in the figure.
Tasks
Task 1
In a rotating reference frame with angular velocity \(\Omega\), state the Bernoulli equation. What are the assumptions used to derive the equation?
Task 2
For the fluid spinning at \(\Omega\), calculate the difference in height of the fluid at points 1 and 2.
Task 3
Calculate the angular velocity of spinning \(\Omega_0\), at which the height of the fluid at point 1 is zero.
Task 4
For the fluid spinning at \(\Omega_0\), calculate the pressure at point 3 in terms of the of the ambient pressure \(p_a\)
5 Mohr’s Circle
The Mohr’s circle is a powerful graphical method used to visualise and analyze the stress state at a single point within a body. Consider two points in a continuum, \(\mathbf{x}_1\) and \(\mathbf{x}_2\), where the stress states are given by
\[ \begin{aligned} \mathbf{S}_1 = \begin{bmatrix} 3 & 1 \\ 1 & 2 \end{bmatrix} \end{aligned} \]
and
\[ \begin{aligned} \mathbf{S}_2 = \begin{bmatrix} 5 & 0 & 0\\ 0 & -1 & -4\\ 0 & -4 & -1 \end{bmatrix} \end{aligned} \]
Tasks
Task 1
We want to derive a non-parametric representation of Mohr’s Circle in 2D.
The transformation rule for second order tensors is given as
\[ \sigma' = \mathbf{R} \, \sigma \, \mathbf{R}^T \]
where \(\mathbf{R}\) is the rotation matrix
\[ \mathbf{R}(\theta) = \begin{pmatrix} cos(\theta) & -sin(\theta) \\ sin(\theta) & cos(\theta) \\ \end{pmatrix} \].
Proceed as follows:
Let \(\sigma_x' = \sigma_n\) be the normal stress and \(\tau_{xy} = \tau_{n}\) be the shear stress. Derive a parametric representation of \(\sigma_n\) and \(\tau_n\) using the transformation rule given above.
Simplify the expression by using trigonometric rules (Hint: we aim for expressions of \(2\theta\))
Obtain a non-parametric expression of a circle, \((x-a)^2 + (y-b)^2 = r^2\).
Task 2
Compute for \(\mathbf{S}_1\) using Mohr’s circle:
The maximum shear stress and the angle \(\theta\)
The maximal normal stress and the angle \(\theta\)
The minimal normal stress and the angle \(\theta\)
Task 3
Compute the principal stresses \((\sigma_1, \sigma_2, \sigma_3)\) of \(\mathbf{S}_2\). We can now construct Mohr’s circle in 3d as follows:
Draw a coordinate system \((\sigma_n, \tau_n)\)
Mark the principal stresses on the abscissa
Draw three circles, each circle passes through two of the three principal stresses marked on the abscissa.
Task 4
Compute for \(\mathbf{S}_2\) using Mohr’s circle:
The maximum shear stress.
The maximal normal stress and the correpsonding directional vector.
The minimal normal stress and the corresponding directional vector.
6 Weak formulation in linear elastics
We consider the momentum balance in static equilibrium
\[ \nabla \cdot \mathbf{\sigma} = \mathbf{f} \]
where \(\mathbf{\sigma}\) is the Cauchy stress tensor and \(\mathbf{f}\) are external forces applied to the system. This equation is the basis for many structural engineering problems.
It is the goal of this exercise to derive a weak formulation of the problem in a simplified setting. We assume Hooke’s law in it’s simplest form and linearly relate the stress \(\mathbf{\sigma}\) to the strain \(\mathbf{\epsilon}= \nabla \mathbf{u}\), where \(\mathbf{\epsilon}\) is the strain, \(\mathbf{u}\) is the displacement and \(c \in \mathbb{R}\) being the stiffness.
\[ \mathbf{\sigma} = c \; \mathbf{\epsilon} = c \; \nabla \mathbf{u} \]
Given the rectangular domain \(\Omega\) shown in the figure with Dirichlet boundaries denoted with \(\Gamma_D\) and Neumann boundaries denoted with \(\Gamma_N\), we state the following strong formulation of our linear elastics problem
\[ \begin{cases} \begin{aligned} & c \nabla \cdot \nabla \mathbf{u} = \mathbf{f} \\ & \mathbf{u} = \mathbf{g} \\ & \nabla \mathbf{u} \cdot \mathbf{n} = \mathbf{h} \\ \end{aligned} \quad \begin{aligned} & \quad \text{on} \; \Omega \\ & \quad \text{on} \; \Gamma_D \\ & \quad \text{on} \; \Gamma_N \end{aligned} \end{cases} \]
Note
Be precise about the dimensions of each object in your derivations. Note that \(\mathbf{u} \in \mathbb{R}^2\) for this two dimensional example.
Tasks
Task 1
Derive the weak formulation of the problem assuming that the displacement \(\mathbf{u} \in H^1(\Omega)\) and test function \(\phi \in H^1(\Omega)\). Make sure to differentiate between the Dirichlet and Neumann boundary term. Why is the implementation of the Dirichlet term in this weak formulation difficult?
Note
Recall the very common notation that a function space with subscript zero, e.g. \(\mathbb{V}_0\) denotes compact support. For our purposes, this means that the function vanishes at the boundary.
Note
Recall from your PDE class that \(H^1\) denotes a Hilbert space.
Tip
Recall from your PDE-class that a weak formulation is derived by multiplying your problem with a test function (e.g. denoted as \(\phi\)) and integrated over the domain \(\Omega\). Often, integration by parts is employed to reduce the demands on the funtional space of the solution.
Your weak formulation always needs to be stated in a way similar to the following:
Find \(\mathbf{u} \in \mathbb{V}\) such that \[ \begin{cases} & \text{YOUR } \\ & \text{PROBLEM} \\ & \text{STATEMENT} \\ \end{cases} \] \(\forall \phi \in \mathbb{W}\)
Where \(\phi \in \mathbb{W}\) denotes the test function we used to derive the weak form.
Task 2
A typical solution strategy to solve the issue with the Dirichlet boundary condition is to ‘lift’ the problem. We will do this step by step
Construct a (simple) function \(\mathbf b(x, y)\) that fulfills \(\mathbf b(0, y) = \mathbf g_0(y)\) and \(\mathbf b(1, y) = \mathbf g_1(y)\).
Define a new function \(\mathbf w = \mathbf u - \mathbf b\)
Consider a suitable function space for \(w\). What changed compared to the function space of \(u\)?
Substituve \(\mathbf w\) into the strong formulation of our problem. Do not forget the boundary conditions.
Write down the weak formulation for the problem statement derived in 4. Make sure to choose an approriate (and useful) function space for the test function as well.
Assume you are given the solution \(\mathbf w^*\) for the weak formulation you posed in 5. Write down the solution for \(\mathbf u\).
7 Stress vector
You are given the following stress vector (tractions)
\[ \begin{pmatrix} \mathbf{t}^{n_1}, & \mathbf{t}^{n_2}, & \mathbf{t}^{n_3} \end{pmatrix} = \begin{pmatrix} \sqrt{2} \begin{pmatrix} 8 \\ 3 \\ -1 \end{pmatrix}, & \sqrt{2} \begin{pmatrix} -2 \\ -3 \\ 5 \end{pmatrix}, & \begin{pmatrix} -3 \\ 2 \\ 0 \end{pmatrix} \end{pmatrix} \]
and the normal vectors
\[ \begin{pmatrix} \mathbf{n}_1 & \mathbf{n}_2 & \mathbf{n}_3 \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 & 0 \\ -1 & 1 & 0 \\ 0 & 0 & \sqrt{2} \end{pmatrix} \]
Tasks
Task 1
Compute the Cauchy stress tensor \(\mathbf{\sigma} \in \mathbb{R}^{3x3}\) by directly solving the linear system associated with \(\mathbf{t}^i = \mathbf{\sigma} \mathbf{n}^i\) for \(i=\{1,2,3\}\).
Note
You can use your computer to solve the system if you want.
Task 2
Now we want to solve for the stress tensor without the need to solve a linear system. You can do so by changing the coordinate system such that
\[ \mathbf{\tilde{t}} = \mathbf{\tilde{\sigma}} \mathbf{n} = \mathbf{\tilde{\sigma}} \mathbf{e} \]
where \(\mathbf{n}\) is an arbitrary vector of unit length and \(\mathbf{e}\) is a cartesian unit vector. Then the rows of \(\mathbf{\tilde{\sigma}}\) are given \(\mathbf{\tilde{t}}\).
8 Divergence Theorem
Recall that the Divergence Theorem for a vector field can be stated as follows: Let \(\Omega \in \mathbb{R}^3\) be a regular domain with boundary \(\partial \Omega\) where \(\mathbf{n}\) is the outward pointing normal on \(\partial \Omega\). Given an arbitrary vector field \(\mathbf{v}\), the relation
\[ \int_{\Omega} \nabla \cdot \mathbf{v} \: d \: \Omega = \int_{\partial \Omega} \mathbf{v} \cdot \mathbf{n} \: d \: A \]
Tasks
Task 1
We are now interessted in gaining a better understanding of the physical meaning of the theorem.
For that matter, proceed as follows:
Write down the Divergence Threorem while explicitly keeping track of the dependence on the spatial coordinate \(\mathbf{x} \in B\), where \(B = (\mathbf{y}, \delta)\) indicates a ball with center point \(y\) and radius \(\delta\).
Now we would like to rewrite the divergence term in an expansion around the center point \(\mathbf{y}\) and higher order terms
Simplify the integral by introducing \(vol(B)\).
Isolate the divergence term and take the limit for \(\delta \rightarrow 0\).
Interpret the final relation in your own words.
Task 2
Given \(\Omega\) and \(\partial \Omega\) as before, we now consider an arbitrary second-order tensor \(\mathbf S \in \mathbb{R}^{3\times3}\). Prove the Tensor Divergence Theorem
\[ \int_{\Omega} \nabla \cdot \mathbf S \: d \: \Omega = \int_{\partial \Omega} \mathbf S \cdot \mathbf n \: d \: A \tag{1}\]
Tip
Multiply (dot product) Equation 1 with an arbitrary (constant) vector. Make use of the Divergence Theorem to proof the Tensor Divergence Theorem.